function EPSP_TC()
% Analyze the time course of EPSP changes
% Analysis includes: amplitude, slope, apparent IPSP, Ihold, and Vrmp.
% 11/1/99. Updates to clean up.
%

%try
sf = getmainselection;
if(sf > 0)
    pflag = getplotflag;
    QueMessage('CC analysis', 1); % clear the que
    for i = 1:length(sf)
        EPSC_TC2(sf(i), pflag, i);
    end;
end;
% catch
%    watchoff;
%    QueMessage('Error in EPSC_TC Analysis routine', 1);
% end;


function EPSC_TC2(sf, plot_flag, i)

global VOLTAGE CURRENT ALLCH
global CONTROL DFILE
% first generate empty output arrays

err = 0;
EPSPresult=[];

QueMessage('EPSP_TC - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on


% this next routine reads all of the data and information we need,
% including calculating some standard stuff that is put into DPAR.
% you will want DPAR.time especially - it is the time base that goes
% with the data.

[DFILE, DPAR, err] = ana_setup(DFILE, sf);
if(err ~= 0)
    return;
end;
isrc = 1;

if(~isempty(ALLCH))
    VOLTAGE = ALLCH{isrc};
    CURRENT = ALLCH{isrc+1};
end;


%-----------------------------Begin analysis--------------------------------

% the analysis consists of 4 parts:
% 1. measure holding current, resting membrane potential, tau, and find any
% spikes in the record (EPSPs in the vicinity of spikes should not be
% counted, and in a perfect LTP experiment, they would not occur while
% measuring epsps or epscs.
%
% 2. We then measure the EPSCs or EPSPs. 
% 3. We do a few more computations and store the results in the database
% 4. We plot the data. 



QueMessage('EPSP_TC - Measuring Ih and rmp');

% measure holding current and "rmp"
% also measure Rin with s1 pulse after ho (ts2)
% and tau by fitting to s1 pulse.
for i = 1:DPAR.records
    hold_cur(i) = mean(CURRENT(i,1:DPAR.tstep(1))); % save in array for later usage
    RMP(i) = mean(VOLTAGE(i,1:DPAR.tstep(1)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(RMP);

if(strcmpi(CONTROL(sf).mode, 'cc'))
    dmode = 0;
    data = VOLTAGE;
else
    dmode = 1;
    data = CURRENT;
end;

% Now smooth the voltage out a bit
QueMessage('EPSP_TC - Smoothing');
for i = 1:DPAR.records
    fsamp = 1000/DPAR.RATES(i); % get sampling frequency
    fco = 1000;		% cutoff frequency in Hz
    wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
    if(wco < 1) % if wco is > 1 then this is not a filter!
        [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
        vsmo(i,:) = DigitalFilt(b, a, data(i,:)')'; % filter all the traces...
    else
        vsmo(i,:) = data(i,:);
    end
end

%vsmo = voltage; % in case you don't want filtering
QueMessage('EPSP_TC - finding spikes');

% find spikes in the records. We can't analyze an EPSP if spikes are present.
if(dmode == 1)
    vsmo = -vsmo;
end;
[rmp_spike, first_isi, nr_spikes, spike_train] = find_spikes2(DFILE, vsmo,  0, sum(DPAR.tsw), DPAR.spike_thresh);


srmp_list = find(isnan(rmp_spike) | rmp_spike >= DPAR.tsw(1)); % NO spikes in rmp window
srmp_nlist = find(~isnan(rmp_spike) | rmp_spike < DPAR.tsw(4)); % spikes IN rmp window

sp_list = cell(2,1); sp_nlist = cell(2,1); s2_list = []; s2_nlist = [];
TM_SPIKE1=[]; TM_SPIKE2=[];

for j = 1:size(DPAR.tpsp, 2)
    for i = 1:DPAR.records
        swin = floor([1:DPAR.tpspi(1, j,i) DPAR.tpspi(2, j,i):DPAR.tpspi(3, j,i)]);
        if(any(vsmo(i,swin)>DPAR.spike_thresh))
            sp_list{j} = [sp_list{j} i];
        else
            sp_nlist{j} = [sp_nlist{j} i];
        end;
    end;
end;
if(dmode == 1)
    vsmo = -vsmo;
end;

TM_SPIKE1=DPAR.TM([sp_nlist{1}]);
TM_SPIKE2=DPAR.TM([sp_nlist{2}]);

QueMessage('EPSP_TC - measuring Rin, Ihold, Tau');
rtdur = floor((DPAR.tstep(2,1)-DPAR.tstep(1,1))*0.8); % duration of window to use for rin measurement
tr1 = DPAR.tstep(1,1) - 2; % baseline for voltage before step (in msec)
if(tr1 < 1)
    tr1 = 1;
end;
tr2 = DPAR.tstep(2,1)-rtdur;  % start of measurement in the pulse
% get baseline voltage and ihold current
vbase = zeros(DPAR.records, 1);
ibase = zeros(DPAR.records, 1);
for i = 1:DPAR.records
    vbase(i)= mean(VOLTAGE(i,1:DPAR.tstep(1,i)));
    ibase(i)= mean(CURRENT(i,1:DPAR.tstep(1,i)));

    if(rtdur > 0) % make sure valid values are available
        dv = mean(VOLTAGE(i, tr1:DPAR.tstep(1,i))) - mean(VOLTAGE(i, tr2:DPAR.tstep(2,i))); % voltage difference
        di = mean(CURRENT(i, tr1:DPAR.tstep(1,i))) - mean(CURRENT(i, tr2:DPAR.tstep(2,i)));
        Rin(i) = 1000*dv/di;

        [a1 a2 taum] = expfit(DPAR.time(i, DPAR.tstep(1,i)+3:DPAR.tstep(2,i))-DPAR.time(i, DPAR.tstep(1,i)), ...
            VOLTAGE(i, DPAR.tstep(1,i)+3:DPAR.tstep(2,i))); % time constant meausurement
        Tau(i) = taum;
    else
        Rin(i) = NaN;
        Tau(i) = NaN;
    end;
end;

Tau(find(Tau <= 0.2)) = NaN; % eliminate too-small tau's
Tau(find(Tau > 5*DPAR.tsw(2))) = NaN; % eliminate too-large taus

% NOW find EPSP amplitudes
if(dmode == 0)
    base = vbase;
    sign = 1; % in cc, measure EPSPs
else
    base = ibase;
    sign = -1 % in vc, measure EPSCs (negative going)
end;


%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

QueMessage('EPSP_TC - measuring EPSP/EPSC');
k = 0;
V_epsp_1 = zeros(DPAR.records, size(DPAR.tpsp, 2));
I_epsp_1 = V_epsp_1;
I_ipsp_1 = V_epsp_1;
T_epsp_1 = V_epsp_1;
T_s_epsp_1 = V_epsp_1;
S_epsp_1 = V_epsp_1;
O_epsp_1 = V_epsp_1;
V_ipsp_1 = V_epsp_1;
S_ipsp_1 = V_epsp_1;
T_ipsp_1 = V_epsp_1;
T_s_ipsp_1 = V_epsp_1;


V = 0;
for j = 1:size(DPAR.tpsp, 2) % for each window to analyze
    for i = sp_nlist{j} % for all the traces WITHOUT spikes detected... 
        pspbase = mean(vsmo(i, 1:DPAR.tpspi(1, j ,i)));
        base(j,i) = pspbase;
        k = k+ 1;
        [V I_epsp_1(k, j)] = max(sign*(vsmo(i,DPAR.tpspi(2, j,i):DPAR.tpspi(3, j,i))'-pspbase));
        V_epsp_1(k,j) = mean(sign*(vsmo(i,(I_epsp_1(k)-2+DPAR.tpspi(2, j,i)):(I_epsp_1(k)+2+DPAR.tpspi(2, j,i)))'-pspbase)); % smooth it out w/ 5 points
        T_epsp_1(k,j) = DPAR.ZT(i); % make matching time base

        [S_epsp_1(k,j) slx] = max(diff(sign*(vsmo(i,DPAR.tpspi(2, j,i):DPAR.tpspi(3, j,i))))/DPAR.RATES(i)); % slope     
        S_epsp_1(k,j) = mean(diff(sign*(vsmo(i,slx-2+DPAR.tpspi(2, j,i):slx+2+DPAR.tpspi(2, j,i))')))/DPAR.RATES(i); % slope - average around the pont
         if(~isempty(DFILE.shkl))
            O_epsp_1(k,j) = (slx+DPAR.tpspi(2, j,i))*DPAR.RATES(i)-DFILE.shkl; % save onset time.
        else
            O_epsp_1(k,j) = (slx+DPAR.tpspi(2, j,i))*DPAR.RATES(i)-DPAR.tstim(1); % save onset time.
        end;
        T_s_epsp_1(k,j) = DPAR.ZT(i); % timing for slopes...

        [V I_ipsp_1(k,j)] = min(sign*(vsmo(i, DPAR.tpspi(3, j,i):DPAR.tpspi(4, j,i))')); % ipsps.
        V_ipsp_1(k,j) = -mean(sign*((vsmo(i,(I_ipsp_1(k,j)-2+DPAR.tpspi(3, j,i)):(I_ipsp_1(k)+2+DPAR.tpspi(4, j,i)))'-pspbase))); % smooth it out w/ 5 points
        T_ipsp_1(k,j) = DPAR.ZT(i);

        [S_ipsp_1(k,j) slx] = max(diff(sign*(vsmo(i, DPAR.tpspi(3, j,i):(I_ipsp_1(k,j)+DPAR.tpspi(4, j,i)))'))); % slope - find point of max slope on RISING phase of IPSP
        S_ipsp_1(k,j) = -mean(diff(sign*(vsmo(i,slx-2+DPAR.tpspi(3, j,i):slx+2+DPAR.tpspi(3, j,i))')))/DPAR.RATES(i); % slope - average around the pont
        T_s_ipsp_1(k,j) = DPAR.ZT(i); % timing for slopes...

    end;
end;

V_epsp_1(find(V_epsp_1<0)) = NaN; % clear out useless information
S_epsp_1(isnan(V_epsp_1)) = NaN;
S_epsp_1(find(S_epsp_1<0)) = NaN;
O_epsp_1(isnan(S_epsp_1)) = NaN;

V_ipsp_1(find(V_ipsp_1<0)) = NaN;
S_ipsp_1(isnan(V_ipsp_1)) = NaN;
S_ipsp_1(find(S_ipsp_1<0)) = NaN;

TMBL=DPAR.ZT(srmp_list);
VBL=vbase(srmp_list)'; % our baselines are also abstracted this way...
IBL=ibase(srmp_list)';
RINL = Rin(srmp_list);
TAUL = Tau(srmp_list);

% determine records for EPSP averaging and measurement, etc.
% basic window is 3 - 0 min prior to first valve change, and last 3 min's
% of the associated shock protocol for that valve change.
% we use the exclusion lists to actually generate a short list to use.

sw1 = []; sw2 = [];
rv1 = []; tv1 = []; iv0 = []; rv0 = []; ep1l = [];
if(~isempty(DPAR.n_valve))
    sw1=DPAR.sw_valve(find(DPAR.n_valve==1)); % records when valve switched to 1.
    sw2=DPAR.sw_valve(find(DPAR.n_valve==2)); % records when valve switched to 2.
end;
if(~isempty(sw1) & ~isempty(sw2))
    ep1 = [sw1(1):sw2(1)-1];
    if(length(sw1) > 1)
        ep2 = [sw2(1):sw1(2)-1];
    else
        ep2 = sw2(1);
    end;

    rv1=sw_valve(find(diff(n_valve)==1)+1); % record and time at which first occured
else
    sw1=0; sw2 = 0; ep1 = 0; ep2=0; rv1=1;
end;
tv1=DPAR.ZT(rv1(1)); % use only the first one...
iv0=find(DPAR.ZT>(tv1-3*60)); % convert from mins to secs, get the time
tv0=DPAR.ZT(iv0(1));
rv0=iv0(1);
if(sw1 == 0 & sw2 == 0)
    sw1 = 1;
    sw2 = max(sp_nlist{1});
end;
ep1l=find(sp_nlist{1} >= sw1(1) & sp_nlist{1} < sw2(1));

k = 0;
for i = sp_nlist{1}(ep1l)
    k = k + 1;
    %	fprintf('[1] i: %d  maxv: %6.1f\n', i, max(vsmo(i,epsp_b1(i):epsp_e1(i))));
    EPSP1(k,:)=vsmo(i, DPAR.tpspi(1,1,i):DPAR.tpspi(4, 1,i)) - base(1,i);
    T_EPSP1(k,:)=DPAR.time(i, DPAR.tpspi(1,1,i):DPAR.tpspi(4, 1,i));
end;
% rotate before plotting;
EPSP1 = EPSP1';
T_EPSP1 = T_EPSP1';
pair = 0;
if(pair)
    ep2l=find(s2_nlist >= sw1(1) & s2_nlist < sw2(1));
    k = 0;
    for i = sp_nlist{2}(ep2l)
        k = k + 1;
        fprintf('[2] i: %d  maxv: %6.1f\n', i, max(vsmo(i,epsp_b2(i):epsp_e2(i))));
        EPSP2(k,:)=vsmo(i, epsp_b2(i):epsp_e2(i)) - base(2,i);
        T_EPSP2(k,:)=DPAR.time(i, epsp_b2(i):epsp_e2(i))-time(i,epsp_b2(i));
    end;
    % rotate before plotting;
    EPSP2 = EPSP2';
    T_EPSP2 = T_EPSP2';
end;

QueMessage('EPSP_TC - analysis complete');

%-----------------------------Prepare for plotting--------------------------------

% for plotting, do baselines/stddev.
bl = find(T_epsp_1 < 5 & ~isnan(V_epsp_1));
mv_e1 = mean(V_epsp_1(bl));
std_e1 = std(V_epsp_1(bl));
pstd_e1 = mv_e1 + std_e1;
nstd_e1 = mv_e1 - std_e1;

ms_e1 = mean(S_epsp_1(bl));
stds_e1 = std(S_epsp_1(bl));
pstds_e1 = ms_e1 + stds_e1;
nstds_e1 = ms_e1 - stds_e1;
tplot = [T_epsp_1(1)+0.1,T_epsp_1(end)-0.1];

mv_i1 = mean(V_ipsp_1(bl));
std_i1 = std(V_ipsp_1(bl));
pstd_i1 = mv_i1 + std_i1;
nstd_i1 = mv_i1 - std_i1;

ms_i1 = mean(S_ipsp_1(bl));
stds_i1 = std(S_ipsp_1(bl));
pstds_i1 = ms_i1 + stds_i1;
nstds_i1 = ms_i1 - stds_i1;

% save the results
if(DFILE.mode >= 5)
    PSP_TC.zt = DFILE.ztime; % we need ztime to align multiple traces...
else
    PSP_TC.zt = DFILE.ztime / (1000);
end;

% store the results.
PSP_TC.T_epsp_1 = T_epsp_1;
PSP_TC.V_epsp_1 = V_epsp_1;
PSP_TC.T_s_epsp_1 = T_s_epsp_1;
PSP_TC.S_epsp_1 = S_epsp_1;
PSP_TC.T_ipsp_1 = T_ipsp_1;
PSP_TC.V_ipsp_1 = V_ipsp_1;
PSP_TC.O_epsp_1 = mean_var(O_epsp_1); % just return the mean max slope time.

% baseline measurements for comparison...
PSP_TC.T = TMBL;
PSP_TC.Vm = VBL;
PSP_TC.Rin = RINL;
PSP_TC.Ih = IBL;
PSP_TC.Tau = TAUL;
PSP_TC.ftime = DFILE.ftime; % save the ftime value too....
% and average EPSC during the last 2 minutes of the sample period.

CONTROL(sf).PSP_TC = PSP_TC;

%----------------------------- plot if figure is set--------------------------------

if plot_flag>=0
    h = findobj('Tag', 'EPSP_TC'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'EPSP_TC', 'Name', 'EPSP Time Course Analysis', 'NumberTitle', 'off');
    end
    figure(h); % otherwise, select it
    clf; % always clear the window...
    fsize = 7;
    msize = 3;
    tmax = max(time);

    %  plot valve on top
    subplot('Position', [0.07, 0.95, 0.6, 0.03]);
    [xs, ys]=stairs(DPAR.TL, DPAR.VL);
    plot(xs, ys, '-k', 'LineWidth', 2);
    set(gca, 'FontSize', fsize);
    ylabel('Valve');
    set(gca, 'XTickLabelMode', 'Manual');

    %plot Max voltage (or current) for EPSP
    if(dmode)
        yax_units = 'pA'
    else
        yax_units = 'mV'
    end;
    subplot('Position',[0.07,0.78,0.6,0.16]);
    plot(T_epsp_1, V_epsp_1, 'ks', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    hold on;
    plot(tplot, [mv_e1 mv_e1], '-k', 'linewidth', 1);
    plot(tplot, [pstd_e1 pstd_e1], '-k', 'linewidth', 1);
    plot(tplot, [nstd_e1 nstd_e1], '-k', 'linewidth', 1);

    if(pair)
        plot(T_epsp_2, V_epsp_2, '-ro', 'Markersize', msize, 'Markerfacecolor', 'r');
    end;
    hold off;
    set(gca, 'FontSize', fsize);
    ylabel(sprintf('EPSP/C Amplitude (%s)', yax_units));
    set(gca, 'XTickLabelMode', 'Manual');
    %grid;
    u=get(gca, 'YLim');


    %	plot Slope of EPSP
    subplot('Position',[0.07,0.61,0.6,0.16]);
    plot(T_s_epsp_1, S_epsp_1, 'k^', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    hold on;
    plot(tplot, [ms_e1 ms_e1], '-k', 'linewidth', 1);
    plot(tplot, [pstds_e1 pstds_e1], '-k', 'linewidth', 1);
    plot(tplot, [nstds_e1 nstds_e1], '-k', 'linewidth', 1);
    if(pair)
        plot(T_s_epsp_2, S_epsp_2, '-ro', 'Markersize', msize, 'Markerfacecolor', 'k') % data in red
    end;
    hold off
    set(gca, 'FontSize', fsize);
    grid;
     ylabel(sprintf('EPSP/C slope (%s/ms)', yax_units));
    set(gca, 'XTickLabelMode', 'Manual');

    %	plot Min (IPSP) voltage
    subplot('Position',[0.07,0.44,0.6,0.16]);
    plot(T_ipsp_1, V_ipsp_1, '-ko', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    set(gca, 'FontSize', fsize);
    grid;
    ylabel(sprintf('IPSP/C Amplitude (%s)', yax_units));
    xlabel('T (min)');
    u=get(gca, 'YLim');

    %	plot Slope of IPSP
    subplot('Position',[0.07,0.27,0.6,0.16]);
    plot(T_s_ipsp_1, S_ipsp_1, 'k^', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    hold on;
    plot(tplot, [ms_i1 ms_i1], '-k', 'linewidth', 1);
    plot(tplot, [pstds_i1 pstds_i1], '-k', 'linewidth', 1);
    plot(tplot, [nstds_i1 nstds_i1], '-k', 'linewidth', 1);
    if(pair)
        %		  plot(T_s_ipsp_2, S_ipsp_2, '-ro', 'Markersize', msize, 'Markerfacecolor', 'k') % data in red
    end;
    hold off
    set(gca, 'FontSize', fsize);
    grid;
     ylabel(sprintf('IPSP/C slope (%s/ms)', yax_units));
    set(gca, 'XTickLabelMode', 'Manual');


    subplot('Position',[0.07,0.22,0.6,0.05]);
    vsp1=15+mean(VBL)+0*TM_SPIKE1; % just set to mean value + 15 mV
    vsp2=15+mean(VBL)+0*TM_SPIKE2; % just set to mean value + 15 mV
    plot(TMBL, VBL, '-ko', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    %  text(TM_SPIKE1, vsp1, '1', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'color', 'blue', 'FontSize', 6);
    text(TM_SPIKE2, vsp2, '2', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'color', 'red', 'FontSize', 6);
    set(gca, 'FontSize', fsize);
    set(gca, 'YLim', [-85 -45]);
    grid;
    ylabel('RMP (mV)');
    set(gca, 'XTickLabelMode', 'Manual');

    subplot('Position',[0.07,0.165,0.6,0.05]);
    plot(TMBL, IBL/1000, '-ko', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    set(gca, 'FontSize', fsize);
    u=get(gca, 'YLim');
    m = max(abs(u));
    set(gca, 'YLim', [-m m]);
    grid;
    ylabel('Ihold (nA)');

    subplot('Position',[0.07,0.11,0.6,0.05]);
    plot(TMBL, RINL, '-ko', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    set(gca, 'FontSize', fsize);
    u=get(gca, 'YLim');
    set(gca, 'YLim', [0 u(2)]);
    grid;
    ylabel('RIN (mohm)');

    subplot('Position',[0.07,0.045,0.6,0.06]);
    plot(TMBL, TAUL, '-ko', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
    set(gca, 'FontSize', fsize);
    u=get(gca, 'YLim');
    set(gca, 'YLim', [0 u(2)]);
    grid;
    ylabel('Tau (msec)');


    % plot the EPSPs analyzed in the top part.

    subplot('Position', [0.72, 0.5, 0.25, 0.5]);
    plot(T_EPSP1, EPSP1, '-k');
    hold on;
    if(pair)
        plot(T_EPSP2, EPSP2, '-r');
    end;
    hold off;

    % plot textual information abstracted from analysis...
    subplot('Position',[0.72,0.0,0.25,0.40])
    axis([0,1,0,1])
    axis('off')
    text(0,0.9,sprintf('%-12s R[%d:%d]\nProtocol:%-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
    text(0,0.8,sprintf('Sol:%-12s  Gain:%4.1f\nLPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 7);

    orient landscape
    drawnow
    % control printing and closing of window for automatic runs
    % f = 1 creates plot and leaves it up
    % f = 2 creates plot but closes it when done
    % f = 3 creates plot and prints it and then closes it
    if (plot_flag > 0)
        print -dljet3;
    end
    if plot_flag == 2
        close
    end
end

